Effect of boundary conditions on diffusion in two-dimensional granular gases 
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We analyze the influence of boundary conditions on numerical simulations of the diffusive prop- 
erties of a two dimensional granular gas. We show in particular that periodic boundary conditions 
introduce unphysical correlations in time which cause the coefficient of diffusion to be strongly de- 
pendent on the system size. On the other hand, in large enough systems with hard walls at the 
boundaries, diffusion is found to be independent of the system size. We compare the results obtained 
in this case with Langevin theory for an elastic gas. Good agreement is found. We then calculate 
the relaxation time and the influence of the mass for a particle of radius R s in a sea of particles 
of radius Rb- As granular gases are dissipative, we also study the influence of an external random 
force on the diffusion process in a forced dissipative system. In particular, we analyze differences in 
the mean square velocity and displacement between the elastic and inelastic cases. 

I. INTRODUCTION 

The interactions among grains and between grains and the boundaries influence profoundly the macroscopic behavior 
of granular systems. To study such complex many body systems, numerical simulations are frequently used where 
one of the most important ingredients are the collision laws introduced to treat interactions^. For dilute assemblies of 
grains one can use molecular dynamics algorithms where periodic boundary conditions are usually used. If the system 
is initially in a square box, a particle going out on the left re-enters the system on the right. We will show below 
that this kind of boundary condition modifies the general dynamics of the grains and introduces large correlations 
in time. This changes the diffusive behavior of the grains. In this paper, we propose an alternative approach to 
calculate numerically the coefficient of diffusion, accurately and with only very small finite size effects. To validate 
our methods, we compare our results for an elastic gas with the Langevin theory. 

As granular gases are dissipative it is necessary to feed energy into the system to keep the particles agitated. To 
thermalize the system, we choose a random acceleration added to each grain at regular time step intervals dt. Our 
final goal in this paper is to study the dependence of the dynamic properties of the granular gas on the mode used to 
force the system. This work is a first step towards understanding the diffusion process in a binary system composed 
of two grain sizes. The system considered here is composed of one particle, s, of radius R s in a sea of particles of 
radius Rb- The particles are spheres constrained to move in a plane and which interact along their equators so that 
the system is two dimensional. The system considered here is dilutc_with a packing fraction of 30 %. The simulations 
are done with the molecular dynamics algorithms (time step driveru and event driveno). 

To characterize the diffusive behavior, we focus on the mean square displacement of the s particle. It is well known 
that for a 2D gas, the integral of the auto-correlation function does not converge!]. This means that the mean square 
displacement does not vary linearly with time. Therefore, strictly speaking, we cannot define a diffusion coefficient in 
2D. However, we show that in a limited range of time, in the stationary state, the mean square displacement can be 
approximated by the linear function: 

< (r(t + t ) - f(t )) 2 >oc ADt (1.1) 
where D can be interpreted as a diffusion coefficient. All quantities are expressed in arbitrary units. 

II. CHOICE OF BOUNDARY CONDITIONS 

In this section we show that periodic boundary conditions introduce strong correlations and therefore alter the 
diffusion process. 

A. Periodic boundary conditions 

Consistent with common practice, we have used periodic boundary conditions to simulate a system of identical 
spheres R s = Rb = 0.5. Initially the particles are placed randomly in a square box of length L. The number of 
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particles is calculated for each system depending on L, R s and the packin g fra ction. Periodic boundary conditions are 
applied in both directions. In this case, for elastic or forced gases (section |TV|), we have observed a strong dependence 
of D (or of the mean square displacement) on the system size. 

In Fig. [I], we have plotted the mean square displacement < (f(t + to) — f(to)) 2 > (Fig. 0a) and J C(t)dt (Fig. |l]b), 
both calculated in the stationary state, as function of t. C(t) is the normalized autocorrelation function: 

rM _ < v(to+t)v(t ) > -<v(t ) >- 

° {t) - < v(t ) 2 > - < Wo) > 2 ■ ( } 

First, we note that the mean square displacement, at large time, varies linearly with time as expected but the slope 
of the curve, i.e. the diffusion coefficient, increases with system size. We show, in the inset to Fig. [l]a, that this 
dependence on L appears already at short time, when < (r 2 (t + t ) — f(t )) 2 L 2 . This feature can be also 

observed in J C(t)dt, which is proportional to the diffusion coefficient. Similarly, we observe that the relaxation 

o 

time r r (i.e. C(r r ) ~ 0) increases with size. In summary, the bigger the system is, the longer the characteristic ti 
t t and the larger the diffusion coefficient D are. We recall that such dependence has been observed by Alder et a 
They proposed the following law for the dependence of D on the number of particles, N, 

D(N) — £)(oo)(l — 2/N). (2.2) 

However their numerical simulations do not support this conjecturei since they fail to observe any saturation of D 
for large systems. In addition, they found strong correlations in the velocity field characterized by the presence of 
vortex flow pattern at the microscopic scale. Our results confirm the lack of convergence for D with system size. In 
addition, this variation of D with L is also observed in the case of inelastic collisions. 

Another important remark is in order. If the system size is, for example, 60 (with 1400 particles of radius R = 0.5), 
the characteristic time r r is found to be around 20 which represents about 200 collisions for a particle. This means 
that a particle needs to undergo 200 collisions to lose completely the memory of its past. According to the Boltzmann 
theory this time should be limited to only a few collisions. Therefore we cannot accept this result as a valid macroscopic 
description of a gas. It is worth noting that the same results are found for both, the time step driven and the event 
driven algorithms. 

We now discuss some points helpful for understanding the problem. Initially, each particle has a random velocity 
drawn from a Maxwellian distribution. We shift the linear and angular momenta so that the system has zero center 
of mass momentum and zero angular momentum relative to the center of mass. We find, however, that, although the 
system keeps its center of mass at rest throughout the simulation, the system is no longer isotropic, its moment of 
inertia becoming that of an ellipsoid. Let I(t) be the inertia matrix of the system. Its two eigenvalues A n and X p are 
related via 

N 

A„ + A p = TO^r2(t), (2.3) 

where the sum is over all N particles each of mass m. Following A„ and X p in time shows that the system takes an 
ellipsoidal form (A„ < A p ). We have found, as well, an anisotropy in the diffusion tensor D(t) defined from I{t) as: 

<"> 

As example we show, in Fig. || the two eigenvalues D\ and D 2 of D as functions of t for a particular periodic system. 
Clearly, Di and D 2 are very different for all t. For all systems we studied, we found two different diffusion coefficients 
which depend strongly the system size. We were not able to find how these values scale with L. 

In addition, we have found that, contrary to its initial condition, the system starts to rotate. This fact is put in 
evidence by calculating the two eigenvectors u n and u p of I(t). These two (perpendicular) vectors rotate in space 
and, most importantly, they keep the same direction of rotation for a long time (~ r r ). We suspect that this rotation 
induces an anomalous temporal correlation of velocities. One should point out that this rotation phenomenon seems 
similar to that observed by Alder et al in their simulations with similar periodic boundary conditions. 

We strongly believe that the use of periodic boundary conditions is responsible for this anomalous correlation. 
These boundary conditions present another inconvenience which is connected with the rotation of the system: The 
square geometry of the system does not permit the conservation of distances between two particles when the system 
is rotating. In Fig. [| we show that after a rotation of 9 the distance dij between particles i and j can be drastically 
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changed by the rotation if one of them goes through the boundary. Because the distance between particles is not 
conserved by rotation, the interaction potential used in the algorithm, which depends only on the relative positions 
of particles dy, is itself not invariant under rotation and so the angular momentum of the system is not conserved. 
Effectively the total angular momentum is fluctuating as one can see in Fig. |]. Every time a particle goes through the 
boundary, its angular momentum, /*, changes sign. Consequently, the change in angular momentum is AL Z = — 2l\. 
AL Z is always proportional to L (the system size) and the total number of particles, N, is proportional to L 2 . However 
it appears that the fluctuations of L z get bigger with system size (see Fig. ^). 

The use of periodic boundary conditions amounts to replicating the system on a square lattice. There are, therefore, 
several identical systems which interact through the boundaries. The rotation observed in our system is then extended 
to all these systems and can create some shear stress, due to frustration of rotation, between neighboring systems. 

These boundary conditions can have other consequences on the dynamics of granular systems. For example, during 
the simulation of a .cooling state the system evolwes towards clustersQ whose orientation depends on the type of 
boundary conditions!! In other similar simulationscl it was shown that there exist large spatial correlations between 
particles where the velocities stay correlated over a distance of about L/2. 



B. Reflecting boundaries 



In the case of reflecting boundaries the system is rotationally invariant leading to better behavior of the mean 
square displacement. However, in this case, the mean square displacement is limited at long time by the system size. 
To circumvent this problem, we proceed as follows. The test particle s is initially put at the center of system at t = 0. 
The evolution of the position and velocity of this test particle are then followed until it reaches the boundary of the 
system in time t w . We then repeat this many times collecting statistics for many test particles with different initial 
velocities. 

The mean square displacement is calculated over 500 such trajectories and limited to time smaller than the smallest 
t w . In this case, as one can see in Fig. ||, there is no dependence of the mean square displacement on the system size. 
Therefore, we can now trust the results of our numerical simulations. 

We recall that the integral of the velocity correlation function does not converge in 2D. However, in a limited range 
of time (see Fig. H), the quantity < (r(t + t Q ) — r(t )) 2 > can be approximated by a straight line and D calculated 
according to Eq. ( |i.l|) . Therefore, the estimate of D with this method is an approximation. 



III. DIFFUSION IN AN ELASTIC GAS 



We first validate our algorithm using reflecting boundaries for an elastic gas (i.e., where the collision between 
particles are elastic). Then, we compare the numerical results with those given by the Langevin equation. Indeed, 
near equilibrium, the dynamics of s can be described approximately, by a Langevin equation: 

C?Vj (t) ,. _ 

_^i = - 7Vi( t ) + r i (t) ) (31) 
<r i (t)r J -(t') >= q 5 id s(t-t'). 

where i denotes the two direction x and y. Integrating Eq. ( |3.l| ), the dependence of the mean square velocity on time 
is simply given by: 

v 2 (t) = v 2 (0)e~ 2 ^ + i(l - e~ 2 ^). (3.2) 

7 

In this paper v denotes the instantaneous velocity of one particle and v 2 the mean square velocity averaged over the 



different s trajectories. We can rewrite Eq. (3.2) using the mean square velocity in the equilibrium state w 2 (oo) = q/j: 



v 2 (t) = v 2 (oo) + {v 2 (0) - w 2 (oo))e- 27t . (3.3) 

where I/7 corresponds to the relaxation time. 

For example, if R s » Rb or, equivalently, m s >> mt (m s .b is the mass of the particle of radius R s ,b), 1/7 is very 
large: the collision of s with a light particle b will not affect strongly the velocity of s. So a great number of collisions 
is needed before s reaches its equilibrium state. Knowing the total kinetic energy of the system Et ot (which is given by 
the initial velocity of each particle), we can easily calculate the square velocity in the equilibrium state ?j 2 (oo) = q/j 
(using the Boltzmann distribution law for elastic gases). Using the simulations to calculate v 2 (t) for the s particles for 
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initial conditions which are very different from the stationary state, i.e. v 2 (0) ^ v 2 (oo) and comparing with Eq. ( |3.3| ), 
we obtain the relaxation tim e for a big (heavy) particle (Fig. [?]). Note that we are looking for agreement near the 



equilibrium state, where Eq. (3.1) is valid. Indeed the dissipation term 7 must depend on both velocities v 2 and v 2 , 
as we will sh ow. 

Equation ( |3.l[ ) also gives the mean square displacement as a function of time, 

< im - mr >= k 2 - q -) il ^ !l + ^ - ^ - (3-4) 

Comparing Eq. (3.4) and Eq. ( fbl] ) at large time, the coefficient of diffusion is seen to be D = - 2 ^°°- ) . In Fig. 0, we 
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compare the theo retic al mean square displacement, Eq. (3.4), with the numerical one, obtained for the case R s = 3R b 



Note that in Eq. (3.4) all the parameters are known. Clearly, the agreement is very good. This confirms that even for 
a large test particle, the motion is well described by the simple Langevin equation. This observation, while reasonable, 
is not trivial since the limited size of our system and the radius of the particles are comparable to the mean free path. 

We can now present a theoretical calculation of 7 which describes dissipation in the Langevin equation for all pairs 
(R s ,Rb). This will allow us to compare the theoretical values with the numerical ones as a function of R s . 

The value of 7 depends on both velocities, v s and Vb- To estimate theoretically the value of 7, we consider the 
deviati on, due to a collision, of the particle s moving at v s in the x direction. The dissipative term —7^ appearing 
in Eq. ([O]), in the x direction, can therefore be formally written as 

v' s .x-v s 

-jv s =< > ui c v s (3.5) 



where v'„ is the velocity after the collision and u> c is the rate of collision. The symbol <> in Eq. (3J5) corresponds to 
the average over all collisions between the s particle and the b ones. To calculate the different terms, we proceed as 
follows. We consider the collision of s with a particle b moving at a velocity Vb- The collision is characterized by two 
angles: the angle between (f* s — rj) and the x axis, and tp the angle between v*b and the x axis. Then, for such a 
collision, illustrated in Fig. 0, we can calculate theoretically v' s (9, ip), the final velocity of the s particle. 
For elastic collisions, the projection of v' s (9, ip) on the x direction is given by, 

v' s (9, (p).x — — v s cos (9) H VbCos(9)cos(9 — ip) + v s sin (9), (3-6) 

m s + mb m s + mb 

with the collision taking place only if 

v s cos(9) — VbCOs(9 — ip) > 0. (3-7) 
Integrating over tp and taking into account Eq. (|3.7j) we can write 
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J £ v (v s cos(9) — VbCos(9 — ip))v' s (9, (p)xdip 

< v> s {9)-x > V = T, ' ( 3 - 8 ) 

/ £ v (v s cos(9) — VbCos{9 — ip)) 







where £ v is the Heavyside function. We found for Eq. ( |3.8| ) two solutions depending on the velocities. If v s < Vb, we 
have 



< v<M.x >„= ^-^ Vs cos 2 { 9) + Vs sin\0) - ^Wglf^ (3 . 9) 
m s + mb (tt — Up) [m s + nib ) 



for all 6 (0 < 8 < tt) and with 9 p — cos 1 (v s cos(9)/vb). For the second case, v s > Vb, there is a critical angle, 
9 C = cos~ 1 (vb/v s ), such that for n — 9 C < 9 < tt + 9 C the collision does not take place. In this case the solution of 
Eq. (O) is 



< v'(6).x >„= — —v s cos 2 (9) + v s sin 2 (9) for < 9 < 9 C , 

m s + m b 

^/„x -, m.-mi 2mbVbCos(9)sin(9 v ) . „ n (o.luj 

< v>(0)£> v = — -v B cos 2 (9) +v s sin 2 (9) - . n w w for6» c < 9 < tt - 9 C . 

m s + m b (tt - 9 p )(m s + m b ) 
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We call v(9) the mean relative loss of velocity, v{6) — <v s ( 6 ^ x> r Vs ) averaging only over the angle (p. In Fig. || we 
show v{&) for the particular case R s = 0.25 and Rb = 0.5, which means that v s > Vb- Note that in our calculation, 
the terms v s and Vb correspond to the averaged values with respect to the appropriate Maxwellian distribution. To 
obtain the mean value, v, of v, we average by integrating numerically over 9. 

To conclude the calculation of the dissipative term, —~fv s , we have to estimate, using Eq. (|3.5|), thc-pollision 
frequency which also depends on the velocities of the two particles. A similar calculation of v can be doneEl In-the 
stationary state where v 2 and v 2 are constant and the distributions of the velocities are Maxwellian one can useO 

uj c = xV^{Rs + R b )d^vl + vl (3.11) 

where d is the density of b particles and \ is a correction factor which corresponds to the local radial distribution 
around the s particle. 

In Fig. |9[ we compare, for different values of R s , the diffusion coefficient found from the simulation with the 
theoretical value, v 2 (oo)/2j, predicted by the Langevin equation combined with our analytical calculation of 7. The 
theoretical calculation of 7 agrees very well with the simulation results. We recall that I/7 corresponds to the 
characteristic time for the diffusive behaviour. It is important to notice that 7 can be approximated by tu c only when 
R s << Rb- Effectively, the calculation for m s ~ gives v = — 1. Larger s particles need to suffer more than one 
collision to lose memory of their previous condition. For m s ~ 00, £> is found equal to zero. Using these methods, 
we find for the elastic monodisperse case (R s = Rb) that relaxation (decorrelation) takes place after about three 
collisions. 

The agreement between numerical results and theoretical predictions allows us to confirm our numerical algorithm. 



IV. FORCED SYSTEM 



In a real granular system dissipation occurs through collisions, a fact that must be taken into account. Experimental 
mechanical properties of grains (restitution and friction coefficients) and collision lawst 2 ! are used in our simulations. 
The collisions between grains and the walls are treated with the same inelastic properties. Due to dissipation, we .need 
to feed energy into the system to maintain the particles agitated. To accomplish this, we choose random heatingE3E-3: 
At every time step St we give a random acceleration, T]i(t), in both spatial directions to each particle. The equation 
of motion can now be written formally as: 

< FP(t)Fj(t') >= m 2 5 tJ S(t - t') v 2 . l4 - ij 

F? is the collision force acting on a particle of mass m. We chose the random acceleration, Ff /m, to be independent 
of the mass of the particle. It is given by a Gaussian noise of variance rfa. 

At long time, the loss of energy due to collisions and the. gain due to F* balance each other such that the system 
reaches a steady state out of equilibrium. It can be showntj that the velocity distribution in this steady state is well 
described by a Maxwellian. 



A. Stationary state 

In the stationary state energy loss and gain balance exactly. The energy loss per unit time, T, for the s particle, 
can be expressed as: 

T = P(m s ,m b )uJ c m s v 2 , (4.2) 

where P{m s , nib) is the relative energy loss of particle s due to collisions. Clearly as for v, T must depend on the mass 
of the particle and on the two velocities v s and Vb- On the other hand, the gain of energy due the stochastic force is 

^m s [v 2 (t + St) - v 2 (t)} = m s r)l5t. (4.3) 

In the steady state of the monodisperse system (R = R s = Rb, and u 2 (oo) =constant), we find, using Eq. ( |3.1l| ), the 
following scaling for v 2 (oo) : 

- 2 M^o) 2/3 (4 . 4) 

V (OO) OC T c 
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We checked these two scaling laws numerically (see Fig. [LC]) and obtained the correct exponent 2/3 for the various 
coefficients of restitution used in the contact laws. We have also verified the predicted dependence on r c for different 
values of R. The good agreement between theory and simulation indicates that we can describe the system by 
macroscopic continuous equations if 8t « r c . As we explain elsewherecJ, the term P{m,m) (in the monodisperse 
case) is independent of mass and velocity, because all particles are identical. This value of Mijn 1 m) was found equal 
approximatively to 0.145 for the mechanical properties corresponding to acetate spheres^ 2 ]. We can then, in the 
case of a mono-disperse system, predict the dependence of the mean square velocity on the various parameters and, 
consequently, characterize the stationary state. For the bi-disperse case, the calculation is more complicated. Indeed 
the loss of energy depends on the two types of colliding .particles and also on the different coefficients of restitution 
and friction introduced in the collision laws. As we showEj the dependence of P(m s ,mb) on v s /vb is not trivial. 

In this paper we limit ourselves to the effect of the thermalization mode (or random force) on the diffusion coefficient. 
To this end, we will compare in the following section the simulation results for D with v 2 (oo)/2j from the Langevin 
equation. 



B. Diffusion of one particle 



To estimate D, we use reflecting boundaries and the same method explained in section [IB . We consider here 
the bi-disperse case (a single particle of radius R s in a sea of particles of radius Rb). As we have not yet found a 
theoretical expression for P(m 8 ,mb) for this case, we use for the mean square velocities the values obtained from the 
simulations which are shown in Fig. [ll|a. Note that ry 2 , has been chosen such that the value of v 2 is the same as in the 
previous section. We see that u 2 (oo) first decreases with R s for R s < Rb but then increases when R s > Rb- Because 
of dissipation and the random acceleration, the repartition of the energy with the mass i s no longer proportional to 



X/m s - In all ca ses i t is possible to calculate the mean collision frequency for s with Eq. (3.11) and the associated 7 
value with Eq. (3.5). We can then calculate the relaxation time r r for all couples (R s ,Rb) used. In Fig. [ll| b we show 



the diffusion coefficient D, and the relaxation time t t -in the inset- as functions of R s . The behavior of v 2 (po) strongly 
modifies the curve of r r and D versus R s compared to the elastic case. Note that the relaxation time represented in 
Fig. [ll] clearly increases as R s increases. 

We have seen in the elastic case that D = |-. In Fig. lib, we show the numerical results for D as a function of R s 



and the corresponding values given by |-. One can see clearly that the external noise modifies the dynamics of the 
granular gas and in particular the diffusion coefficient, D. The numerical value of D is found to be larger than that 
obtained by the corresponding random walk. Indeed, at short time, due to the random force, v 2 (t) is not constant. 
Between two collisions v 2 (t) increases linearly with t. Starting with the equation of motion of particle s (between two 
collisions) , 

^ = *(*), («) 
and with the initial conditions Xi(0) and Vi(0), we can calculate mean square displacement 



< (ssi(t) - x,(0)) 2 >= ( / dh Vi (0) + / 77i(ti)dti \ J dt 2 I vi(0) + / m (t 2 )dt' 2 } ) . (4.6) 

\i=0 \ / t 2 =0 

In two dimensions, this yields for the interval between two collisions 

< (r(t) - r(0)) 2 ) >= v 2 (0)t 2 + ^t 3 . (4.7) 

On the other hand, in the case of a random walk (or elastic collisions) the mean square displacement at short time 
scales as t 2 . This difference explains the disagreement between D and " 3°°) _ As the velocity changes between two 
collisions the probability of collision is increasing with time too. The calculation of the coefficient of diffu sion is not 



easy in this case, due to the correlation between the velocity and the probability of collision (see Eq. p.ll|) . For very 
small particles, if one approaches relaxation by the time of a new collision, i.e. v ~ — 1, this calculation should be 
possible. Indeed we can assume that the velocities before and after a collision are not correlated and have the same 
distribution (molecular chaos). We .can then compute the mean square displacement, knowing the dependence of the 
collision probability on the velocityEx With this assumption we improve the estimate of D for the smallest R s . But 
for bigger particles we have seen that the velocities stay correlated over many collisions and we can no longer use 
molecular chaos. 
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V. CONCLUSIONS 



We have presented here some general results about the diffusion process in an agitated granular gas. We first 
showed that the boundary conditions used in the simulations are of crucial importance. Indeed, periodic boundary 
conditions introduce artificially strong temporal correlations which alter the macroscopic properties of the gas. If we 
ensure that no correlations are induced by the algorithm, for example by using reflecting boundaries, the numerical 
results obtained for an elastic gas can be described very well by a Langevin equation. We have presented a theoretical 
calculation of the relaxation time which allows us to predict the diffusion coefficient in all cases studied. This was 
not a priori intuitive since the radius of the particles is of the order of the mean free path. Finally we have analyzed 
the influence of uniform heating (a random acceleration) on dissipative gases. We have shown that heating influences 
the dynamics at short time. This is evident through the value of the diffusion coefficient which is different from that 
expected from the Langevin description. We are now applying with success these results to the diffusion process in a 
granular mixture consisting of two type of grains (differing by mass or size) in equal proportion. 
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FIG. 1. Dependence of the mean square displacement on the system size, (a): < (r(t + to) — r{to) > as a function of t. From 
bottom to top the system size is 20, 40, 30, 50, 60. (b): Integral of C(t) as a function of t for the same system. From bottom 
to top the system size is 20, 40, 30, 50, 60, respectively 



FIG. 2. Di and D2 vs time for a typical monodisperse case with periodic boundaries. 



FIG. 3. illustration of the non-conservation of the distance between two particles. 
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FIG. 4. Total angular momentum L z (t) vs t for two system lengths and the same particle density. (Full line): L = 20, 
(Dashed line): L = 50. 

FIG. 5. < (r(t) - r(0)) 2 > as function of t using reflecting boundaries. Superposed (as in figure are the results for system 
sizes: 20, 30, 40, 50, 60. 



FIG. 6. A particle of mass m a colliding with a particle of mass nib: definition of the angles 9 and ip. 
FIG. 7. Comparison between numerical results and Langevin approximation. R s = 1.5 and R p = 0.5. (a) dependence of 



the mean squared velocity on t; Q : numerical result, full line: Fit using e 7 accor ding to Eq. (3.2). (b) the mean squared 



displacement, : numerical result, full line: Theoretical prediction according to Eq. (3.3) using for 7 the value obtained from 
Fig. §a. 

FIG. 8. v{0), the mean relative loss of energy per collision in the 9 direction, for R s = 0.25, Rb = 0.5 and v% = 25 

FIG. 9. Coefficient of diffusion for different values of R a . Rb = 0.5. : Numerical values obtained by simulation. Full line: 
theoretical values calculated from - i 00 - 1 



FIG. 10. (a) v 2 (oo > versus t]q for different coefficients of restitution. : e n = 0.87, e s = 0.4, fj, — 0.25; A: e n — 0.4, 
e s = 0.4, fi = 0.25. (b) u 2 (oo) versus the mean time between collisions t c . 



FIG. 11. (a) v 2 (oo) vs R s (Rb = 0.5). (b) Coefficient of diffusion D as a function of R s . Q: Numerical values obtained 

^(00 

w 



from simulation. Full line: Corresponding values given by v . The insert shows r r versus R 
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